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Abstract. Using the methods of general relativity Lindquist derived the radiative transfer equation 
that is correct to all orders in v/c. Mihalas developed a method of solution for the important 
case of monotonic velocity fields with spherically symmetry. We have developed the generalized 
atmosphere code PHOENIX, which in 1-D has used the framework of Mihalas to solve the radiative 
transfer equation (RTE) in 1-D moving flows. We describe our recent work including 3-D radiation 
transfer in PHOENIX and particularly including moving flows exactly using a novel affine method. 
We briefly discuss quantitative spectroscopy in supernovae. 
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INTRODUCTION 

In the 1970s Dimitri Mihalas and collaborators developed the equations of radiative 
transfer and numerical techniques for their solution [1, 2, 3, 4, 5, 6] in relativistically 
expanding atmospheres that are required for the quantitative spectroscopic modeling 
of supernovae (SNe). Over the last twenty years our group has developed PHOENIX 
to include most of the physics that is needed to calculate supernova atmospheres. While 
many of the techniques in PHOENIX are modem [7, 8, 9, 10, 1 1], the underlying solution 
of the equations follows the method developed by Mihalas [1]. While this method works 
well for spherically symmetric flows, it is difficult to extend it to the case of fully 3-D 
atmospheres even with the restriction of monotonic velocity fields (see the contribution 
by J. I. Castor, this volume). Here we describe a method of using an affine parameter 
in order to calculate the solution along geodesies (straight lines in flat spacetime). This 
method is straightforward, exact to all orders in v/c, and can be generalized to arbitrary 
flows and include the effects of curved spacetime. 



MOTIVATION 

Why study supernovae? Supernovae are common events that occur about once a second 
out in the volume out to redshift z ~ 1. Supernovae play a major role in galactic 
nucleosynthesis. Supernovae inject energy into ISM, and trigger star formation. And 
finally, SNe la make good standardizable candles. It is important to emphasize while 
the standard candle relation is purely empirical [12, 13] it is understood theoretically. 
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at least qualitatively as being due to higher temperatures due to nickel mass variations 
[14, 15, 16] and opacity variations [14, 17]. 

Observationally, supemovae are classified by their spectra and there are a large num- 
ber of classifications which are described in the review article by Filippenko [18]. The- 
oretically, there are two supernova mechanisms, core collapse and thermonuclear. Core 
collapse occurs at the end of the life of a massive (M > 8 M©) star, when the iron core 
collapses to nuclear matter density and then bounces. Exactly how the shock gets out 
of the iron core is a subject of current research. The wide variety in observed spectra 
is thought to be due to variations in the progenitor (does it have an intact hydrogen en- 
velope or has it lost some or all of its envelope either in a wind or via interaction with 
a companion). Thermonuclear supernova refers to the thermonuclear burning of a C-i-O 
white dwarf that accretes mass from a companion until it reaches the Chandrasekhar 
mass, these objects form the class SNe la. Since most, in not all, SNe la explode at the 
same mass this naturally accounts for the relatively small amount of variation in their 
intrinsic brightness at maximum light. Since SNe la are quite bright they can be found 
far away. The small variation in brightness combined with their high intrinsic brightness 
makes them excellent cosmological probes. 

The small intrinsic variation of SNe la at maximum light can be corrected for using 
the light curve shape method. This follows from the realization that the luminosity 
at peak is correlated with the rate at which the light curve declines from maximum 
light [12]. Refining the original suggestion of Phillips has led to improved light curve 
shape parametrizations [13, 19, 20]. Using SNe la as correctable candles, two groups 
discovered the dark energy [21, 22]. As noted above the light curve shape relations are 
purely empirical and determined using nearby supemovae. Thus, if the distant sample is 
significantly different from the nearby sample this could result in large systematic errors 
in the results for the nature of the dark energy (the existence of the dark energy seems to 
be on pretty solid ground). One way to search for differences is to compare spectra of the 
nearby sample to those in the distant sample. Riess et al. [21, see their Figure 1 1] found 
that the nearby sample was pretty similar to the distant sample when one accounted for 
the variation in the signal to noise. Recently, a modest change with redshift has been 
observed [23]. These authors speculate that this variation is simply due to the enhanced 
star formation rate at high redshifts and thus one is looking at a sample with younger 
progenitors. 

Regardless of whether or not the diversity can be corrected for using purely empirical 
methods, an understanding of the diversity is needed. To first order the diversity in the 
peak brightness has been understood due to a variation in the total amount of nickel 
that is produced in the explosion [14]. This is understood physically as increasing the 
temperature which then leads to the observed variations in the spectra [16]. Figure 1 
shows that simply varying the model temperature does a good job of reproducing the 
observed variation. The spectral sequence also led to the discovery that there are spectral 
indicators that are also well correlated with the brightness at maximum light [16, 24]. 
These spectral indicators can be used as complementary to light curve shape methods 
and since they have relatively small wavelength baselines they are rather insensitive to 
properties of dust in the parent galaxy. Figure 2 shows the definition of the ratio "^fisis 
which Bongard et al. [25] showed can be used by IDEM. Recently, the SNfactory [26] 
showed that a spectral ratio can be defined that reduces the variation in the Hubble 
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FIGURE 1. The observed variation in spectra with brightness at maximum light is well reproduced by 
variation the model temperature in synthetic spectra calculations. 
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FIGURE 2. The definition of the specti-al ratio 9^, 
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diagram to 12%. 



3-D Explosions 

The explosion of a SN la is inherently three dimensional since the white dwarf be- 
comes convective shortly before it ignites and the ignition is likely to be off-center. The 
flame also quickly becomes wrinkled and thus there will be non-spherical compositions. 
Overall, SNe la are basically round, but explosion asymmetries lead to compositional 
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and ionization inhomogeneities. Core-collapse supemovae almost certainly come from 
asymmetric engine. How important the asymmetry is to spectra formation depends on 
how much of the hydrogen or helium envelope has been lost and how late one observes 
the spectrum. Gamma ray bursts are beamed and highly relativistic. And of course AGN 
are asymmetric and may require general relativity. 

DETAILED QUANTITATIVE SPECTROSCOPY 

PHOENIX 

PHOENIX is a generalized model atmosphere code, that solves the generalized stel- 
lar atmosphere problem in static or moving flows using the fully special relativistic 
approach. For supernovae with characteristic velocities of 10,000 km s~^ and maxi- 
mum velocities of up to 60,000 km s~^ the exact special relativistic formulation is pre- 
ferred. PHOENIX is well calibrated on many astrophysical objects including the sun 
[27], cool stars [28], hot stars [29], planets [30], novae [31], and all types of supemovae 
[32, 33, 34, 35, 36]. 

Over the last four years or so we have been developing a fully 3-D version of 
PHOENIX [37, 38, 39, 40, 41, 42]. We have built this up slowly and carefully mak- 
ing sure to do careful code verification along the way. Fig. 3 shows the results from 
a spherically symmetric static test case with scattering. This test was especially useful 
because it showed that our full characteristics method was far superior to a short char- 
acteristics method in terms of reproducing the spherical symmetry. Since this model is 
a sphere in a box, the sphere is surrounded by vacuum and the interpolations needed for 
the short characteristics method did a very poor job of reproducing the correct results 
across the opacity discontinuity. 

Figure 4 shows a toy solar model using periodic boundary conditions. The tempera- 
ture and density structure was taken from the model of Caffau et al. [43] but the opacity 
was taken to be proportional to the density and independent of temperature. The ther- 
malization parameter in this calculation was taken to be e = 1. 

Our general method is described in detail in Refs [37, 38, 39, 40] and is modeled 
on the methods in our 1-D code [9, 10, 11]. The mean intensity J is obtained from the 
source function S by a formal solution of the RTF which is symbolically written using 
the A-operator A as 

J = hS. (1) 

The source function is given by 5 = (1 — e)J+eB, where e denotes the thermal coupling 
parameter and B is Planck's function. 

The A-iteration method, i.e. to solve Eq. 1 by a fixed-point iteration scheme of the 
form 

•Aiew = •'^'^old) "^new = (1 ~ c)'^new + (2) 

fails in the case of large optical depths and small e. The idea of the ALI or operator 

splitting (OS) method is to reduce the eigenvalues of the amplification matrix in the 
iteration scheme [44] by introducing an approximate A-operator (ALO) A* and to split 
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FIGURE 3. Sphere in a box, mean intensity / is plotted for a a slice through the x-y plane for the case 
of a scattering continuum with e — 10^^. 

A according to 

A = A* + (A-A*) (3) 

and rewrite Eq. 2 as 

7„ew = A*S„ew + (A- A*)Sold. (4) 

This relation can be written as [45] 

[1 - A* ( 1 - £)] 7„ew = 7fs - A* ( 1 - e)7oid, (5) 

where 7fs = ASom and 7oid is the current estimate of the mean intensity 7. Equation 5 is 
solved to get the new values of 7 which is then used to compute the new source function 
for the next iteration cycle. 

Matrix computations are required in order to solve Eq. 5. In choosing matrix methods 
one must remember that the inverse of a banded matrix is full. Band matrix solvers 
require large amounts of memory. This is also the case for parallelized band matrix 
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FIGURE 4. The top left panel shows the continuum, the top right panel the line center, the bottom left 
the line wing, and the bottom right a composite image. 
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FIGURE 5. Convergence rates of the 3D transfer for line transfer with plane-parallel test structures 
(label 'PP') and the 3D hydro structure (label 'hydro'). For comparison, the convergence of the A iteration 
for plane-parallel continuum transfer is also shown. 



solvers. Iterative methods (Jordan and Gauss-Seidel) work well and use little memory. 
Ng acceleration is still very useful. Figure 5 shows that the ALO works well and the Ng 
acceleration is evident. 
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MOVING ATMOSPHERES 



For the moving case one must derive the equation of radiative transfer using the tech- 
niques of general relativity. The pioneering work in this field was done by Lindquist [46]. 
We sketch the development of the solution of the relativistic radiative transfer equation 
beginning with the work of Mihalas [1] and describe the work we have done [41, 42]. 



Mihalas' Method 

In his important paper Dimitri Mihalas [1] explained why we want to work in the 
co-moving frame (italics as in original): 

The emissivity rj and opacity x depend upon the angle as well as frequency in 
the inertial frame because of Doppler shifts, aberration, and advection induced 
by the motion of the material in the frame. 

The goal of this section is to rewrite equation (2.1) with all material and 
radiation-field quantities measured in the comoving frame; in that frame both 
the opacity and emissivity are isotropic, and can be related directly to proper 
variables that specify the thermodynamic state of the material. Furthermore, in 
that frame both the scattering properties of the material and the rate equations 
describing the mechanisms populating and depopulating its internal energy 
states are most easily defined. . . 

In our analysis we shall, however, leave both the space and time variables in 
the inertial frame, as this is the only frame in which synchronism of clocks can 
be effected, and further this choice obviates the need to develop a metric for 
accelerated fluid frames [47] which in general can only be done approximately. 
With this choice of frame we can write exact Lorentz transformations for all 
the material and radiation-field quantities and use these to develop a transfer 
equation that will remain valid for relativistic flow in the limit as v/c — > 1. . . 

In Mihalas [1] the characteristic equations are written as coupled ordinary differential 
equations in both the inertial frame spatial coordinate, r, and the co-moving frame 
momentum space coordinate, jU, i.e.. 



4^ = r(M+i8), 

dSM 
dSM 



= 7(1 -m') 



r dr 



From these results of Mihalas a generation of workers in the field have been educated 
that the characteristics are curved in phase space [see Figure 1 of Reference 1]. 
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FIGURE 6. Characteristics are straight Hnes in flat spacetime and they begin at a boundary point and 
traverse the computational volume with a constant direction measured in the observer's frame. 

The Affine Method 

Generalizing the method of Mihalas [1] to 3 spatial dimensions (and thus 6 phase 
space dimensions) is very difficult. Following the methodology of Mihalas [1] we would 
keep all 3 momentum space variables co-moving. Then, in order to keep track of the 
evolution of those co-moving phase space variables along the characteristics (e.g., /i, 
0), we would have relied on the so-called tetrad formalism [48, 49, 50] which is very 
complicated for spacetimes with no symmetry or arbitrary flows. For an example, a boost 
followed by an arbitrary 50(3) rotation is still a valid Lorentz transformation and thus it 
is difficult to define a consistent frame in the case of arbitrary flows. In Chen et al. [42] 
we realized that photons travel along null geodesies (straight lines in flat spacetime), 
which depends only on the background spacetime, not the velocity distribution of the 
flow, and consequently is either known analytically or can be numerically solved before 
we solve the RTF. The so-called characteristics in phase space are just the 'unique' 
lifting of the geodesies in the 4-D spacetime manifold to the 8-D tangent bundle of the 
spacetime manifold [46, 51]. The apparent curvature of the characteristics in [1] is due 
to the selection of phase space coordinates. To avoid the complications induced by the 
use of the tetrad formalism, we found that we could work in a very special co-moving 
frame where only the wavelength of the photon is measured by a co-moving observer 
(which differs from the inertial frame value by only a Doppler factor and is thus easily 
calculated), the other two momentum directions are measured in the observer's frame 
(which in flat spacetime are constants). This is not a giant leap, since in the case of 
Mihalas [1] the transfer equation is solved in a mixed frame, where the spatial variables 
are Fulerian, but the two momentum variables ji and A are measured by a co-moving 
observer. Figure 6 illustrates the characteristics for a Cartesian grid. 
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The above considerations lead to a radiative transfer equation given by 




+ a{s)X 



(6) 



where the notation is the same as in Baron et al. [41]. 

The evaluation of in this mixed frame makes it quite easy to solve the scattering 
problem in the co-moving frame. To get the co-moving , we need to integrate the (co- 
moving) specific intensity (r, n) over the co-moving solid angle element dQ.. Since 
Ix (r, n) is given in terms of the inerdal frame direction n, it is desirable to transform 
the integral into one over the inertial frame solid angle element d£lQ^ this causes no 
difficulty since 



where (r, n) is expressed in the "funny frame" that comes from solving Eq. 6. 

We have implemented this method in spherical coordinates for the case of homologous 
flows [41]. It is important to verify that the code gives the correct results in a case that 
can be tested. Luckily, we can compare the results for a spherically symmetric test with 
those of our well-tested 1-D code, which uses the Mihalas' method and not the affine 
method. Figure 7 compares the results for the mean intensity 7 in the co-moving frame 
of the 1-D code (dots) and the 3-D (solid lines) at each voxel at a given radius for 
differing maximum velocities. No scaling is performed between the two calculations. 
The agreement is excellent and the 3-D code does a very good job of reproducing 
spherical symmetry. This also shows that there is no problem in the interpretation of 
our frame. Figure 8 shows similar results for a scattering line in the co-moving frame. 
Figure 9 shows the flux in the observer's frame reproduces the expected P-Cygni profile. 
Figure 10 shows that the code parallelizes quite well in terms of resolution. The scaling 
test has each processor working on the same number of characteristics as the number of 
momentum space angles is increased. The 14% cost is reasonable given the increased 
communication overhead as one goes from 256 to 16,384 processors. 



PHOENIX now includes tested 3-D fully relativistic radiative transfer with homologous 
flows. Including arbitrary flows is a computational, not an algorithmic challenge [52, 
53] and we are working on this. The next step is to go beyond test problems to a 
full production code. We expect progress from the inter-comparison of new datasets 
of nearby supemovae with 3-D hydro models and synthetic spectroscopy. The affine 
method is extremely easy to implement and can be applied to radiation hydro codes. 
Happy birthday Dimitri! 



da = (7[l-j8-n])-2jno 

= fis)-^dao. 



(7) 



Using Eq. 7, in the co-moving frame can be expressed as 




(8) 



SUMIVIARY 
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FIGURE 7. The results of 1-D calculations are compared with 3-D calculations for j3, 
(0.03,0.33,0.67,0.87). 
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FIGURE 8. The mean intensity of a line in the co-moving frame for a line with e =0.1, Pmax = 0.03. 
The solid lines are the 3-D result and the plus signs are the 1-D results. 
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FIGURE 9. The line shown in Fig 8 transformed to the observer's frame gives the expected P-Cygni 
profile. 
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